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HIGHLIGHTS 


•  We  model  the  thermal  and  electrical  dynamics  of  thermoelectric  power  generating  systems. 

•  Both  transient  and  steady-state  conditions  are  considered. 

•  We  develop  a  computer  program  for  transient  simulations  of  thermoelectric  systems. 

•  The  program  simulates  the  electro-thermal  coupled  effects  that  occur  during  changes  in  the  operating  conditions. 

•  Comparison  of  experimental  and  simulation  results  shows  great  accuracy. 


ARTICLE  INFO 


ABSTRACT 


Article  history: 

Received  9  October  2013 

Received  in  revised  form  14  December  2013 

Accepted  18  December  2013 

Available  online  14  January  2014 


Keywords: 

Thermoelectric 

TEG 

Heat  transfer 
Computer  simulation 
Peltier 
Seebeck 


Recent  interest  in  the  use  of  thermoelectric  generators  (TEGs)  to  recover  waste  heat  in  large-scale  appli¬ 
cations  calls  for  precise  simulation  to  appropriately  design  complicated  and  dynamic  systems.  The  aim  of 
this  work  is  to  develop  a  computer  tool  to  accurately  simulate  the  thermal  and  electrical  dynamics  of  a 
real  thermoelectric  (TE)  power  generating  system. 

The  computer-aided  model  presented  here  is  able  to  accurately  simulate  the  non-linear  electro-ther¬ 
mal  coupled  effects  which  occur  during  changes  in  the  operating  conditions,  e.g.  temperature  or  load 
changes. 

Simulation  results  are  compared  to  experimental  data  obtained  from  a  real  TE  system.  The  comparison 
shows  great  accuracy  both  during  transients  and  in  the  steady-state,  thus  validating  the  model  as  a  reli¬ 
able  tool  to  simulate  TE  generating  systems. 

©  2013  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Thermoelectric  (TE)  devices  can  directly  convert  thermal  energy 
into  electrical  energy  and  vice  versa.  Thermoelectricity,  based  on 
the  fact  that  charge  carriers  can  be  set  in  motion  by  a  difference 
of  temperature,  generally  refers  to  two  main  physical  phenomena: 
the  Seebeck  effect  and  the  Peltier  effect.  The  former  states  that  a 
certain  open-circuit  voltage  is  created  in  a  material  kept  between 
two  different  temperatures.  The  Seebeck  coefficient  a  (V/K)  is  a 
material  property  that  relates  the  open-circuit  thermoelectric  po¬ 
tential  V0c  (V)  with  the  temperature  difference  AT  (K),  or 

Voc  =  ocAT  (1) 

The  Peltier  effect  states  that  a  direct  current  1(A)  passing 
through  a  circuit  of  dissimilar  materials  pumps  thermal  power 
from  one  material  to  the  other: 

Pp  =  nl  =  otlTj  (2) 
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where  n  (V)  is  the  Peltier  coefficient  and  the  last  equivalence  comes 
from  the  Kelvin  relationship;  Tj  is  the  junction  temperature. 

Other  important  phenomena  occurring  in  thermoelectric  de¬ 
vices  are  the  well-known  Joule  heating  and  the  Thomson  effect. 
The  latter  states  that  there  is  reversible  absorption  or  liberation 
of  heat  (in  excess  of  the  joule  dissipation  I2R )  in  a  homogeneous 
material  simultaneously  exposed  to  temperature  gradient  and 
electric  current: 

PT  =  tIAT  (3) 

where  t  (V/K)  is  the  Thomson  coefficient,  defined  as 
rp  dec 

t  —  Tavg  (4) 

The  Thomson  effect  is  usually  much  smaller  than  the  Joule  heat¬ 
ing  [1,2]  and  ts  contribute  can  be  significant  under  large  tempera¬ 
ture  differences  [3,4];  also,  the  Thomson  coefficient  is  difficult  to 
obtain  experimentally,  therefore  it  is  often  neglected  in  literature. 
Its  effect  is  not  included  in  this  article. 

Doped  semiconductors  prove  to  be  the  materials  with  the 
best  thermoelectric  properties.  Multiple  pellets  of  p-  and 
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n-semiconductor  are  connected  electrically  in  series  and  thermally 
in  parallel  to  achieve  and/or  sustain  higher  voltages,  thus  forming  a 
TE  module.  This  paper  will  focus  on  TE  generators  (TEGs);  Fig.  1 
shows  a  3D  model  of  a  TEG  module  where  thermal  energy  is  ap¬ 
plied  on  the  bottom  (“hot”  side).  The  pellets  are  electrically  ser¬ 
ies-connected  by  solder  and  an  Aluminium-based  ceramic  layer 
serves  as  electrical  isolation  and  mechanical  substrate.  The  result¬ 
ing  module  is  quite  robust  and  reliable,  and  operates  without  any 
vibration  or  noise.  It  is  commercially  produced  in  a  wide  range  of 
sizes  from  a  few  millimetres  to  several  centimetres  on  a  side.  Mul¬ 
tiple  modules  can  be  electrically  connected  in  series  or  parallel  in 
order  to  achieve  higher  output  voltages  and  currents.  Heat  is  con¬ 
ducted  through  the  module  while  additional  heat  is  generated 
(Joule)  inside  the  module  and  pumped  (Peltier)  through  the  hot 
and  cold  sides  as  described  by  the  coloured  arrows.  In  TE  power 
generation  the  Peltier  effect  is  parasitic  because  it  reduces  the  tem¬ 
perature  difference  across  the  device,  thus  increasing  the  module’s 
effective  thermal  conductivity.  Thermoelectric  heat  pumps  (cool¬ 
ers  or  heaters)  have  been  used  for  many  years  in  applications  rang¬ 
ing  from  IC  microcoolers,  to  refrigerators  [5],  to  power  stations  [6], 
and  they  are  often  referred  to  as  Peltier  devices.  Due  to  their  low 
efficiency,  often  less  than  5%,  thermoelectric  generators  have  on 
the  contrary  been  used  predominantly  in  military  and  aerospace 
projects  or  for  applications  in  which  cost  is  not  as  important  as 
the  ability  to  reliably  generate  power  in  hostile  or  maintenance- 
free  environments.  However,  interest  in  TEGs  has  recently  in¬ 
creased  because  of  concerns  about  climate  change,  coupled  with 
increasing  performance  and  lower  module  cost.  TEGs  can  effec¬ 
tively  lend  themselves  to  sustainable  applications  of  waste  heat 
recovery  in  which  thermal  energy  is  rejected  to  ambient  and  is 
effectively  free,  e.g.  in  vehicles  [7-10]  and  stoves  [11,12].  TEGs 
can  also  be  used  to  harvest  geothermal  energy  [13]  or  combined 
to  PV,  solar  thermal  or  thermophoto  voltaic  systems  14-16].  More¬ 
over,  recent  advances  in  TE  materials  and  ‘mass-production’  vol¬ 
umes  [17,18]  will  continue  to  lead  to  a  further  improvement  of 
TEGs’  efficiency  and  reduction  of  their  cost,  respectively. 

It  is  of  fundamental  importance  to  carefully  design  large-scale 
systems  in  which  materials  cost  is  of  great  relevance.  However 
TE  systems  are  usually  composed  of  heat  masses,  TEGs  and  power 
and  control  electronics  and  they  are  influenced  by  several  thermal 
and  electronic  phenomena  whose  interaction  is  complex.  More¬ 
over  TEGs  are  often  employed  in  dynamic  environments  which  fre¬ 
quently  undergo  thermal  transients.  Actual  CAD  tools  do  not  yet 
include  the  ability  to  model  thermoelectric  effects  therefore  they 


cannot  be  successfully  used  to  accurately  simulate  the  electro¬ 
thermal  coupled  effects  which  take  place  during  changes  in  the 
system  operating  conditions,  e.g.,  temperature,  power  or  load 
changes.  In  literature  some  research  has  focused  on  this  issue; 
Lineykin  and  Ben-Yaakov  [19]  have  divided  the  TE  module  into  a 
grid  of  spatially  discretized  thermal  circuits  then  transformed  into 
their  electrical  analogies.  The  transient  term  is  included  in  the 
model  as  a  parallel  electrical  capacitor  to  take  into  account  an  addi¬ 
tional  time-related  term  due  to  the  change  in  stored  heat  energy.  A 
similar  approach  has  been  followed  by  Chen  et  al.  [20].  The  accu¬ 
racy  of  the  simulation  is  related  to  the  number  of  cells  used  and 
the  parameters  are  difficult  to  obtain  from  manufacturers’  data 
sheets,  but  most  importantly  they  do  not  offer  a  theoretical  solu¬ 
tion  to  the  problem  because  their  governing  equations  are  based 
on  steady-state  solutions.  A  more  appropriate  approach  would  be 
to  study  the  transient  physical  equations  which  describe  TE  de¬ 
vices  behaviour;  among  these  the  most  important  is  the  heat  equa¬ 
tion.  Al-Nimr  et  al.  [21]  have  already  explored  this  possibility  but 
they  used  fixed  temperatures  as  boundary  conditions  at  the  two 
sides  of  the  TE  module,  assuming  that  those  temperatures  are 
not  varying,  i.e.,  supposing  thermal  isolation  (or  steady-state). 
However,  during  thermal  transients  there  is  exchange  of  heat 
through  the  sides. 

Very  recently  two  very  interesting  works  by  Cheng  and  Huang 
[22]  and  Meng  et  al.  [23]  proposed  two  models  for  thermoelectric 
coolers.  The  former  slightly  overestimates  the  temperature  differ¬ 
ence  in  steady-state,  while  the  second  has  a  maximum  error  of 
4.5  K  over  a  temperature  difference  of  around  37  K  (with  a  current 
input  of  1  A). 

We  have  already  provided  a  mathematical  solution  of  the  one¬ 
dimensional  heat  conduction  equation  for  TE  devices  that  includes 
internal  Joule  heat  generation  and  dynamic  exchanges  of  heat 
through  the  hot  and  cold  sides  in  [24]. 

The  aim  of  this  work  is  to  couple  the  aforementioned  solution 
with  the  other  thermal  and  electrical  phenomena  occurring  in  real 
TE  systems.  The  resulting  physical  model,  described  in  Section  2, 
takes  into  account  the  dynamic  relations  between  the  several  ther¬ 
mal  masses  and  the  most  important  thermoelectric  phenomena 
occurring  in  a  generalised  TE  system. 

The  physical  model  is  then  used  as  the  basic  structure  to  devel¬ 
op  a  computer  tool,  described  in  Section  3,  capable  of  accurate  sim¬ 
ulations  of  the  thermal  and  electrical  dynamics  of  a  physical  TE 
power  generating  system.  The  model  is  created  in  Simulink  and 
Matlab  and  a  comparison  between  experimental  and  simulated 
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Fig.  1.  3D  model  of  a  thermoelectric  generator  showing  the  main  physical  effects. 
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results  is  presented  in  Section  4  to  demonstrate  the  effectiveness 
and  accuracy  of  the  proposed  simulation  model. 

2.  Generic  model  of  a  thermoelectric  system 

Generally  speaking,  a  TE  system  is  composed  of  several  thermal 
masses  which  store  and  exchange  heat  through  conduction  and/or 
convection.  Heat  is  partially  converted  to  or  from  electricity  inside 
the  TE  device  and/or  in  other  elements,  e.g.,  electrical  heaters.  It  is 
very  important  to  optimally  design  the  thermal  interconnection  of 
the  TE  device  with  the  rest  of  the  system  to  guarantee  optimal  per¬ 
formance  [25].  This  paper  ultimately  aims  at  developing  a  model 
that  can  be  used  to  aid  TE  system  design. 

Fig.  2  shows  the  architecture  of  a  generic  TE  power  generating 
system.  The  TEG  is  usually  in  contact  with  a  thermal  mass  on  both 
the  hot  and  cold  side.  Power  is  provided  to  or  removed  from  the 
thermal  masses  in  the  form  of  electrical  or  heat  power,  leading 
to  changes  in  the  thermal  energy  stored  inside  the  thermal  masses. 
Through  conduction  or  convection  part  of  this  energy  is  transferred 
to  and  from  the  TEG  module.  The  TEG  is  modelled  considering  the 
sides  of  the  TEG  separately  from  the  inner  part.  The  heat  equation 
(HE)  deals  with  both  the  heat  conduction  and  generation  (Joule 
heating)  inside  the  TEG.  Additional  heat  is  pumped  at  the  sides 
where  there  is  a  junction  of  two  dissimilar  materials  (Peltier 
effect). 

Part  of  the  energy  flowing  through  the  TEG  is  converted  into 
electrical  power  and  the  process  is  closely  related  to  the  thermo¬ 
electric  effects  described  by  the  Eqs.  (1)  and  (2). 

Conduction  and  convection  equations  provide  the  rate  of  ther¬ 
mal  energy  flowing  from  one  medium  to  the  other  in  time.  The 
thermal  conduction  equation  (Fourier’s  law)  is  written  in  time- 
derivative  form  as: 

dEcond  dT 

~dt r~~KAdx  (5) 

where  Econd  0)  is  the  thermal  energy,  k  (W/mK)  is  the  thermal  con¬ 
duction  coefficient,  A  (m2)  is  the  surface  of  heat  exchange  and  x  (m) 
is  the  direction  of  heat  transfer,  perpendicular  to  A.  The  heat  con¬ 
vection  equation  is 

=  hAAT(t)  (6) 

where  h(W/m2  K)  is  the  heat  convection  coefficient  and 
A T(t)  =  T(t)  —  is  time-dependent. 

Any  lumped  thermal  mass  changes  its  temperature  when 
receiving  or  losing  thermal  energy: 

E  therm  —  mCAT(t')  (7) 


where  mC  (J/K)  is  the  heat  capacity  of  the  thermal  mass,  which  de¬ 
fines  its  ability  to  keep  heat  energy. 

As  indicated  in  Fig.  2,  T^h  is  used  in  this  paper  for  the  temper¬ 
ature  of  the  hot  thermal  mass  and  TH  is  the  temperature  on  the  hot 
side  of  the  TEG,  while  T^c  and  Tc  are  the  temperatures  of  the  cold 
thermal  mass  and  TEG  cold  side  respectively. 

3.  Computer  model 

The  solution  to  the  heat  equation  (HE)  presented  in  [24]  was 
programmed  in  Matlab,  therefore  the  Mathworks  suite  was  chosen 
to  develop  the  simulation  program.  The  model  includes  several 
blocks  and  performs  time-step  simulations.  Simulink  is  an  excel¬ 
lent  choice  of  environment  for  this  task. 

Section  3.2  presents  the  heat  equation  block,  Section  3.3  deals 
with  the  thermal  part  and  3.4  explains  the  electrical  part  of  the 
system. 

3.2.  System  architecture 

Fig.  3  shows  the  architecture  of  the  Simulink  model  developed 
in  this  work  which  will  be  described  shortly. 

The  HE  computes  the  transient  solution  starting  from  an  initial 
condition  of  Tc  and  TH,  and  with  constant  values  for  TocH,  TooC  and 
the  load  current  //oad .  It  can  provide  the  temperature  distribution 
inside  the  TE  device  at  any  instant  in  time,  even  if  for  this  model, 
at  system  level,  we  are  interested  in  the  temperatures  only  at  the 
hot  and  cold  sides  of  the  TEG  module.  The  HE  allows  the  side  tem¬ 
peratures  Th  and  Tc  to  change  depending  on  the  temperatures  T^h 
and  Tooc  of  the  thermal  masses  in  contact  with  the  TEG  through 
conduction  or  convection  media.  The  temperature  difference 
across  these  media  changes  dynamically,  thanks  to  the  Newton’s 
Law  of  Cooling,  as  already  described  in  [24]. 

In  the  solution  of  the  HE  the  temperatures  of  the  thermal 
masses  Tooh,c  are  considered  constant;  however,  in  a  real  system 
during  a  transient  they  vary  depending  on  the  energy  added  or  re¬ 
moved  by  Pin  and  Pout •  Similarly,  the  temperatures  at  the  sides  of 
the  TEG  are  influenced  also  by  the  actual  load  current  which  in 
turn  determines  the  value  of  the  Peltier  effect  (see  Eq.  2).  This  is 
why  in  Fig.  3  the  ‘HE  Block’  is  connected  to  both  to  the  ‘Electrical 
Block’  and  to  the  ‘Thermal  Block’.  The  thermal  block  requires  a 
and  I ioad  form  the  electrical  block  to  deal  with  Peltier  effect.  The 
transient  evolution  of  the  temperature  difference  across  the  TEG 
device  and  the  system’s  dynamic  exchange  of  heat  due  to  Pin  and 
Pout  continuously  influence  the  temperatures  in  the  whole  system. 

Therefore  in  order  to  couple  the  continuous-time  transient 
solution  provided  by  the  HE  to  the  temperatures’  evolution  in 
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Fig.  2.  Architecture  of  a  thermoelectric  power  generating  system. 
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Fig.  3.  Block  architecture  of  the  presented  model.  The  inputs  to  the  different  blocks 
are  written  inside  the  arrows. 


the  remaining  parts  of  the  system  (thermal  masses,  input/output 
thermal  energy  and  electrical  load),  the  model  is  designed  in  the 
discrete-time  domain.  Within  each  time  step  of  duration  tstep  (s) 
the  HE  computes  the  dynamics  inside  the  TEG  for  a  transient  of 
duration  tstep.  The  temperatures  in  the  remaining  parts  of  the  sys¬ 
tem  change  as  described  by  Eqs.  (5)-(7),  in  which  the  derivative 
term  is  substituted  by  slope  terms.  This  direct  examination  in  dis¬ 
crete  terms  aids  the  comprehension  of  the  underlying  physical  ef¬ 
fects  in  the  system. 

To  understand  how  the  model  is  structured,  let’s  consider  two 
successive  time  steps,  corresponding  at  two  program’s  discrete  iter¬ 
ations,  called  i  and  i  +  1 .  It  must  be  noted  that  the  calculations  made 
during  each  iteration  are  computed  directly  using  the  equations 
provided  in  this  paper,  without  any  trial-and-error  or  iterative 
method.  During  the  iteration  i  the  heat  equation  uses  the  electrical 
parameters  and  the  heat  masses’  temperatures  Tooh©  and  T^q*), 
considered  constant  during  tstep,  and  the  hot  and  cold  side  temper¬ 
atures  Tm  and  TC(i)  as  initial  values;  it  computes  the  values  of 
and  rC(i+i)  corresponding  to  their  evolution  after  the  transient  time 
tstep  and  outputs  them  through  a  memory  block  of  duration  tstep,  so 
that  they  are  actually  available  only  at  the  start  of  iteration  i  +  1. 
The  HE  block  is  described  in  Section  3.2.  At  iteration  i  the  model  also 
updates  the  temperatures  of  the  thermal  masses  depending  on  the 
quantity  of  heat  energy  provided  by  Pin  and  Pout  during  tstep,  so  that 
they  are  ready  for  iteration  i  +  1.  The  equations  used  to  update 
these  temperatures  are  described  in  the  Section  3.3. 


3.2.  The  TEG  heat  equation  block 

The  solution  to  the  HE  provided  in  [24]  treats  the  TEG  module  as 
a  whole  block,  i.e.  it  does  not  deal  with  the  different  materials 
(ceramic,  semiconductor,  solder)  separately;  it  also  does  not  divide 
the  semiconductor  materials  into  single  pellets  or  thermocouples. 
This  is  done  by  design  because  the  parameters  for  the  model  are 
obtained  only  from  the  physical  dimensions  of  the  module  and 
from  its  direct  use  in  the  system,  i.e.  by  representing  a  real  system. 

The  HE  code  is  programmed  into  a  Matlab  embedded  equation 
block.  The  input  variables  to  the  block  are: 


•  TEG  geometrical  parameters:  surface  A  and  thickness  L. 

•  Conduction  and  convection  coefficients:  open-circuit  thermal 
conductivity  k  and  thermal  diffusivity  coefficient 


thermal  conductivity  k 
volumetric  heat  capacity  pCp 


of  the  TEG;  thermal  transfer  coeffi¬ 


cients  through  the  hot  and  cold  mediums  pH  and  pc. 


•  Electrical  parameters  of  the  TEG:  internal  resistance  Rint  and 
current  J,oad. 

•  Temperatures  at  the  beginning  of  the  iteration:  T^h,  T^c,  Th  and 

Tc. 

•  Time  step  duration  tstep. 

pc  and  pH  may  vary  if  the  materials  in  contact  with  the  cold  and 
hot  sides  change  properties,  e.g.  increased  convection,  e  can  vary 
with  temperature,  but  in  this  model  it  is  considered  constant. 
The  time  step  size  tstep  is  constant  and  chosen  based  on  the  knowl¬ 
edge  of  the  TE  system;  it  is  a  compromise  between  the  thermal 
time  constant  of  the  transient  evolution  inside  the  TEG  and  the 
temperature  change  in  the  thermal  masses  of  the  system.  Values 
between  1  s  and  3  s  are  usually  chosen  to  satisfactorily  model 
the  physical  system.  For  systems  with  a  very  high  energy  flux, 
e.g.,  large  engine  exhaust  gas  systems,  values  of  under  1  s  might 
be  more  appropriate. 

TooHjooHjcjHjioa^Rint  are  fed  iteratively  at  the  beginning  of 
every  iteration  from  the  rest  of  the  Simulink  model. 

As  already  explained  in  Section  3.1,  the  HE  block  produces  two 
outputs,  Tc  and  THt  that  are  the  temperatures  on  the  cold  and  hot 
sides  after  that  particular  iteration  cycle  of  duration  tstep.  Tc  and  TH 
are  then  passed  through  a  memory  block  so  that  they  are  effec¬ 
tively  available  only  after  tstep,  be.,  at  the  beginning  of  the  succes¬ 
sive  iteration.  The  block  could  be  set  to  calculate  temperatures 
inside  the  TEG  too,  but  these  are  not  of  interest  to  the  aim  of  the 
simulation  (nor,  usually,  to  the  system  designer). 

If  there  are  TEG  modules  connected  in  series  or  parallel  some  of 
the  inputs  required  by  the  HE  will  be  scaled  accordingly  depending 
on  the  total  number  of  modules  Ntot :  Atot  =  NtotA ,  while  Rint  and  //oad 
will  be  provided  by  the  electrical  block,  as  described  in  Section  3.4. 


3.3.  The  thermal  block 


The  most  important  tasks  executed  by  the  hot  side  block  are  the 
update  of  the  temperature  T^h  of  the  hot  thermal  mass  and  the  up¬ 
date  of  the  TEG  hot  side  temperature  TH,  accounting  for  the  Peltier 
effect. 

Consider  now  the  input  power  transferred  to  the  hot  thermal 
mass.  This  can  be  provided  by  an  electrical  heater  or  by  thermal 
conduction/convection.  For  the  electrical  case  it  is  easily  measured 
while  in  the  thermal  case  it  can  be  calculated  using  the  well- 
known  laws  of  heat  conduction  or  convection.  As  shown  in  Fig.  3, 
part  of  the  input  power  goes  into  changing  the  quantity  of  heat  en¬ 
ergy  stored  in  the  hot  thermal  mass,  while  the  remaining  thermal 
energy  is  transferred  (through  the  hot  thermal  mass)  to  the  TEG. 
Heat  losses  can  be  accounted  for  in  another  term  PiossH .  It  can  be 
written: 


D  _  WhChAT^h 

1  in  ~  } 

Lstep 


+  fa  A  AT H  +  PiossH 


(8) 


where  mHCH  is  the  hot  thermal  mass  heat  capacity  J/K,  A T^h  is  the 
change  in  temperature  of  the  hot  thermal  mass  during  tstep,  hH  is  the 
convection  coefficient  (or  for  conduction  the  conduction  coefficient 
over  the  thickness)  of  the  medium  between  the  hot  thermal  mass 
and  the  TEG,  A  is  the  area  of  the  TEG  and  A TH  is  the  temperature  dif¬ 
ference  across  the  medium. 

From  Eq.  (8)  it  would  be  possible  to  obtain  the  new  value  for 
Toon,  however,  we  first  need  to  account  for  the  Peltier  effect  which 
acts  on  the  second  term  on  the  right  of  Eq.  (8).  The  Peltier  effect 
pumps  additional  power  from  the  hot  side  to  the  cold  side,  thus 
effectively  reducing  the  hot  side  temperature  TH  provided  by  the 
HE.  If  we  write  that  PconvH  =  hHAATH  =  faAfT^  -  TH)f  where  TH 
is  the  result  of  the  HE,  then  we  can  correct  the  value  of  TH  by  cal¬ 
culating  its  ‘real’  value: 
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■y  y  PconvH  PpeltH 

1  Hreal  =  1  ooH - .  A - 

where  PPeitH  =  ocITH  from  Eq.  (2).  Now  that  the  ‘real’  value  for  TH  is 
obtained,  from  Eq.  (8)  it  is  possible  to  calculate  T^h  for  the  next 
iteration: 


»H(i+ 1) 


=  TooH(i)  + 


]Pm  P lossH  hp[A  (T ooH(i)  T Hreal ) ]  f step 

hIhCh 


(10) 


It  can  be  noted  that,  depending  on  the  polarity  of  the  second  term 
on  the  right  of  Eq.  10,  TooW(I-+i)  can  be  greater  or  smaller  than 
Toon®-  T Hreal  and  ^^(i+i)  will  be  used  in  the  next  iteration  (i  +  1)  of 
the  HE  as  the  initial  value  of  the  hot  side  temperature  and  as  the 
temperature  of  the  hot  thermal  mass,  respectively. 

For  the  cold  part  of  the  system  similar  considerations  hold  true. 
The  power  flowing  to  the  cold  thermal  mass  is  the  sum  of  the  ther¬ 
mal  power  flowing  through  the  TEG  and  the  thermal  power 
pumped  to  the  cold  junction  by  the  Peltier  effect.  In  a  similar 
way  as  for  the  hot  side  we  update  Tc  to  its  ‘real’  value: 


Tcreal  —  T ooC  + 


hcA(Tc  —  Tqqc)  +  odTc 

hrA 


(11) 


It  is  now  possible  to  obtain  TooC(i+ 1)  as 


1  ooC(i+l) 


=  1 


ooC(i)  ' 


\hc  A  (Tcreal  Too  C(i))  Pout  f^ossc]  Ttep 

_ 


(12) 


where  Pout  is  the  power  removed  from  the  cold  thermal  mass  and 
PCout  takes  into  account  any  thermal  power  lost  to  ambient  on  the 
cold  side. 

Both  the  values  calculated  in  Eq.  (11)  and  (12)  are  passed  to  the 
HE  block  for  the  next  iteration  (i  +  1). 


3.4.  The  electrical  block 


A  TEG  can  be  modelled  as  a  voltage  source  V0c  in  series  with  its 
internal  resistance  Rint ,  even  during  transients,  because  their  elec¬ 
trical  dynamic  response  is  in  order  of  nanoseconds  [26].  Therefore, 
when  the  TEG  is  connected  to  a  load  Rioad  it  can  be  written  that  the 
load  voltage  is 

V load  —  V oc  ~  P-inthoad  (13) 

It  can  be  noted  that  both  the  internal  resistance  and  the  open-cir¬ 
cuit  voltage  can  be  approximated  to  linear  functions  of  the  temper¬ 
ature  difference  AT  across  the  device  [27],  therefore  Eq.  (13)  can  be 
updated  to 

Vioad  =  mVocAT  -  (mRjntAT  +  qRint)hoad  (14) 

As  stated  in  the  Introduction,  this  work  neglects  the  Thomson  ef¬ 
fect;  therefore,  mVoc  is  considered  constant  and  equal  to  the  Seebeck 
coefficient  a  (V/K).  mRint  (Q/K)  and  qR.m  (Q)  are  constant  coefficients, 
too. 

In  order  to  use  Eq.  14  in  the  model,  the  appropriate  values  for 
the  coefficients  mVoc ,  mR.nt  and  qR.m  should  be  obtained  through  an 
electrical  characterisation  of  the  TEGs  used  in  TE  system  [28]. 

The  electrical  block  also  deals  with  series  and  parallel  connec¬ 
tion  of  TEG  modules.  It  is  assumed  that  all  the  modules  are  identi¬ 
cal  and  all  are  subjected  to  equal  thermal  conditions.  If  there  are  Ns 
modules  in  series  and  NP  modules  in  parallel,  for  a  total  of  Ntot,  the 
total  open-circuit  voltage  is  V0ctot  =  NsVoc  and  the  load  current  J/oad 
is  NP  times  the  current  in  a  single  branch.  It  can  then  be  written, 
from  Eq.  (13),  that 

V,oad  =NSVoc-’^J^  (15) 

The  dimension-less  figure  of  merit  ZT  and  the  actual  efficiency  rj  are 
also  calculated  14,9]: 


L  oc2  Tpj  +  Tc 
RintA  k  2 


n  = 


Vloadhoad 
Pin  ~  P lossH 


(16) 

(17) 


In  Eq.  16  the  electrical  conductivity  was  written  as  a  =  L/(RintA)  and 
the  average  temperature  of  the  device  as  TAVG  =  (TH  +  Tc)/2. 

The  electrical  part  of  the  Simulink  model  computes  the  values 
of  Rint  and  Vioad  depending  on  the  actual  temperature  gradient  AT 
and  on  the  electrical  parameters  J/oad,  Ns,  NP.  The  Peltier  term  is  cal¬ 
culated  as  N tot ocTjI Singief  where  Isingie  is  the  current  passing  through  a 
single  TEM.  The  load  current  (or  voltage)  in  a  real  TE  system  is  set 
either  by  a  constant  Rioad  connected  to  the  TEG,  or  by  an  interfacing 
power  electronic  converter  with  Maximum  Power  Point  Tracking 
(MPPT)  [29].  For  flexibility  of  simulation,  the  desired  load  is  passed 
to  the  model  as  an  input,  but  the  model  provides  also  the  maxi¬ 
mum  theoretical  power,  i.e.  when  Vioad  =  Voc/ 2,  to  compare  actual 
electrical  power  output  with  the  maximum  that  can  be  extracted. 


4.  Experimental  and  simulation  results 

In  order  to  test  the  ability  of  the  proposed  Simulink-Matlab 
model  to  simulate  real  TE  generating  systems,  a  transient  experi¬ 
ment  has  been  performed  in  the  laboratory  and  then  simulated 
by  computer.  The  results  for  both  and  a  discussion  about  their 
comparison  is  presented  in  this  Section. 

4.1.  Experimental  test  rig 

The  TE  system  used  in  the  experiment  is  identical  to  that  pre¬ 
sented  in  [28].  The  cold  thermal  mass  is  formed  of  a  water-cooled 
copper  block  and  the  hot  thermal  mass  is  a  copper  block  containing 
one  high-power  high-temperature  electrical  heater  powered  by  a 
DC  power  supply  unit.  The  TEG  output  is  connected  to  an  electronic 
load.  The  TEG  module  is  compressed  between  the  two  blocks  by  a 
mechanical  fixture  which  does  not  introduce  thermal  shocks  or 
short-circuit  the  thermal  path.  The  mechanical  compression  force 
on  the  system  is  sensed  by  a  load  cell  and  thermocouples  measure 
the  temperatures  of  the  water,  the  heater,  and  the  faces  of  the  TEG 
module  (directly  touching  them).  All  the  sensors  are  connected  to  a 
datalogger.  Each  of  the  electronic  instruments  used  is  connected  by 
the  IEEE-488  interface  (GPIB)  to  a  computer  and  controlled  by  a 
program  created  in  the  Agilent  VEE  Pro  software. 

The  copper  blocks  have  a  thermal  capacity  of  approximately 
565  J/I<  (Copper  volumetric  heat  capacity  CVcu  =  3.45  J/cm3  K), 
the  thermal  transfer  coefficient  for  between  the  TEG  and  the  hot 
and  cold  thermal  masses  are  respectively  9100  and  1 1,800  W / 
m2  K  (Copper  thermal  conductivity  kCu  =400W/mK;  the  heater 
is  2.2  cm  from  the  TEG  hot  side,  the  water  pipes  are  1.7  cm  from 
the  TEG  cold  side). 

The  input  power  is  provided  by  electrical  heaters  and  thermal 
energy  is  removed  from  the  cold  side  through  forced  liquid  flow. 
The  mass  flow  rate  of  the  water/glycol  mixture  is  measured  as 
17  g/s  and  the  coolant  mixture  heat  capacity  is  3.89  J/gK.  The  con¬ 
vection  coefficient  was  determined  empirically  and  is  1400  W/ 
m2  K.  Heat  losses  from  the  hot  block  (which  is  insulated  with 
glass-fibre)  are  included  in  the  simulation  as  described  in  [28]. 
The  time  step  size  used  for  the  simulation  is  3  s.  This  choice  is 
based  on  the  fact  that  the  test  rig  is  able  to  vary  the  temperature 
difference  at  a  maximum  rate  of  0.2  K/s,  therefore  a  time  step  of 
3  s  ensures  that  the  temperatures  of  the  thermal  masses  do  not 
change  much  and  at  the  same  time  guarantees  a  good  evolution 
of  the  transient  solution  provided  by  the  heat  equation. 

The  TEG  module  used  is  manufactured  by  European 
Thermodynamics  Ltd.  (model  code:  GM250-1 27-1 4-10)  and  has 
an  open-circuit  thermal  conductivity  k=1.62W/mK,  surface 
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A  =  1600  mm2,  thickness  L  =  3  mm,  thermal  diffusivity 

e  =  le-6  m2/s. 

Electrical  characterisation  of  the  module  used  was  performed  at 
two  different  temperature  differences,  AT  —  50  °C  and 
AT  =  150  °C,  both  with  a  mechanical  force  of  2000  N  on  the  mod¬ 
ule,  corresponding  to  a  pressure  of  1.2  MPa.  Results  obtained  are 
shown  is  Fig.  4.  These  data  are  then  used  to  obtain  the  coefficients 
of  Eq.  (14):  mVoc  =  0.048;  mRint  =  0.0045;  qR.m  =  1.52. 

4.2.  Comparison  and  discussion  of  results 

The  designed  transient  lasts  for  55.5  min  and  is  intended  to  pro¬ 
duce  several  power  and  load  steps.  The  VEE  program  records  all  the 
electrical  data,  temperatures,  mechanical  pressure  and  time.  The 
transient  parameters  are  as  follows,  where  Pin  is  the  DC  electrical 
power  to  the  heater  and  J/oad  is  the  current  value  set  in  the  elec¬ 
tronic  load: 

1.  Pin  =  50  W ;Iioad  =  0  A  (open-circuit)  @  time  =  0  s. 

2.  Pin  =  50  Wfioad  =  0.2, 0.4, 0.6, 0.8, 1  A  @  time  =  1195,  1210, 
1229,  1247,  1286  s. 

3.  =  150  W ;h0ad  =  1  A  @  time  =  1286  s. 

4.  P^  »150W;//Ofld*  1.2, 1.4, 1.6, 1.4, 1.2, 1.0, 0.8, 0.6  A  @  time  = 
2761,  2780,  2795,  2828,  2844,  2877,  2946,  3109  s. 


Fig.  4.  Electrical  characteristic  of  the  TEG  module  used  in  the  experiment,  for  two 
different  temperature  gradients  AT  =  50,  150  °C,  clamped  at  2  1<N/1.2  MPa. 


This  experiment  is  designed  to  test  the  simulation  model  under 
thermal  and  electrical  transients  of  difference  magnitude.  Different 
temperature  variation  rates  were  applied  both  in  open-circuit  and 
at-load,  with  several  electrical  load  steps.  The  Simulink  program 
runs  without  interruptions,  therefore  testing  it  under  different 
thermal  and  electrical  conditions  ensures  that  there  are  no  offsets 
that  would  make  the  results  valid  for  just  one  of  those  operating 
conditions. 

Fig.  5  shows  both  the  experimental  and  simulation  results  in 
the  same  graph,  for  comparison.  In  particular  the  temperature  dif¬ 
ference  and  the  output  electrical  power  are  plotted  versus  the  tran¬ 
sient  time.  During  the  transient  the  temperature  difference  across 
the  TEG  device  varies  from  5  °C  to  almost  150  °C,  leading  to  some 
thermal  expansion  effects  in  the  whole  system.  The  mechanical 
pressure  on  the  system  could  not  be  controlled  during  the  test, 
but  active  control  of  the  clamping  force  will  be  incorporated  into 
a  future  test  rig  to  compensate  for  thermo-mechanical  effects.  In 
previous  work  [28]  it  was  noted  that  different  pressure  loads 
slightly  influence  the  contact  resistances  in  the  system,  thus  affect¬ 
ing  the  thermal  conductivity  and  the  power  produced  by  the  TEG 
module.  Their  effect  should  not  account  for  more  than  a  5%  change 
in  the  results. 

As  it  can  be  appreciated  from  Fig.  5,  there  is  no  delay  in  the  re¬ 
sponse  of  the  simulated  results  with  respect  to  the  experimental 
data.  The  simulation  results  show  very  good  agreement  with  the 
experimental  results,  accurately  tracking  the  electro-thermal  cou¬ 
pled  effects  occurring  in  the  TE  system.  The  Root  Mean  Square  Er¬ 
ror  (RMSE)  is  3.58  °C  for  the  temperature  difference  and  0.2  W  for 
the  power  output.  The  RMSEs  normalised  to  the  values  range 
(NRMSE)  become  2.76%  and  4.55%  for  the  temperature  difference 
and  the  output  power,  respectively.  There  are  no  works  in  the  lit¬ 
erature  that  proposed  a  simulation  model  for  TEG  systems  demon¬ 
strating  similar  prediction  performance  under  varying  thermal  and 
electrical  conditions  as  the  one  presented  in  this  article.  As  a  con¬ 
sequence  it  is  difficult  to  compare  the  data  presented  here  to  other 
simulation  results  reported  in  the  literature. 

It  is  very  interesting  to  note  the  close  effect  that  the  load  cur¬ 
rent  has  on  the  temperature  gradient  across  the  device.  Every  in¬ 
crease  in  Iioad  leads  to  a  decrease  in  temperature  due  to  the 
greater  value  of  the  Peltier  term  of  Eq.  (9).  Conversely,  the  last 
500  s  of  the  test  show  that  a  smaller  load  current  increases  the  rate 
of  change  of  temperature.  The  effective  thermal  conductivity  of  the 
TEG  module  varies  markedly,  depending  on  the  current  flowing 
through  the  module.  MPPT  converters  should  take  this  thermal 


Fig.  5.  Experimental  and  simulated  transient  results  showing  the  temperature  difference  and  the  output  electrical  power  versus  the  time. 
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variation  into  account  when  setting  the  optimum  operating  point 
for  maximum  efficiency  or  power  generation. 

Future  work  will  develop  a  model  of  the  system  that  incorpo¬ 
rates  the  non-linear  effects  introduces  by  a  MPPT  converter.  The 
model  will  be  adapted  to  a  typical  waste  heat  energy  recovery  sys¬ 
tem  from  exhaust  gas  for  automotive  applications. 

5.  Conclusion 

This  paper  presents  an  innovative  and  powerful  computer  tool 
to  accurately  simulate  real  thermoelectric  power  generating  sys¬ 
tems  even  during  transient  conditions.  The  proposed  simulation 
program,  designed  in  Simulink-Matlab,  deals  with  all  the  most 
important  thermoelectric  effects  and  is  able  to  cope  with  both 
thermal  and  electrical  dynamics.  Consequently  it  can  greatly  help 
the  design  phase  of  large-scale  and/or  complicated  thermoelectric 
systems. 

A  comparison  of  experimental  and  simulation  results  shows  the 
accuracy  and  capability  of  the  model,  showing  that  it  can  be  em¬ 
ployed  to  study  transients  and  steady-state  operation  of  real  TE 
systems  with  confidence. 
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